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We develop from first principles Markovian master equations suited for studying the time evolution of a 
system evolving adiabatically while coupled weakly to a thermal bath. We derive two sets of equations in 
the adiabatic limit, one using the rotating wave (secular) approximation that results in a master equation in 
Lindblad form, the other without the rotating wave approximation but not in Lindblad form. The two equations 
make markedly different predictions depending on whether or not the Lamb shift is included. Our analysis 
keeps track of the various time- and energy-scales associated with the various approximations we make, and 
thus allows for a systematic inclusion of higher order corrections, in particular beyond the adiabatic limit. We 
use our formalism to study the evolution of an Ising spin chain in a transverse field and coupled to a thermal 
bosonic bath, for which we identify four distinct evolution phases. While we do not expect this to be a generic 
feature, in one of these phases dissipation acts to increase the fidelity of the system state relative to the adiabatic 
ground state. 



I. INTRODUCTION 



Recent developments in quantum information processing, in particular theoretical HI and experimental [2] proposals for 
adiabatic quantum computation, have generated considerable renewed interest in the old topic of quantum adiabatic dynamics 
fl3]|4l. While much work has been done on rigorous formulations of adiabatic approximations for closed quantum systems, e.g., 
Refs. l5l- [T2l . adiabatic evolution in open quantum systems is still a relatively unexplored topic. In this regard master equations 
governing the evolution of a quantum system with a time-dependent Hamiltonian coupled to an external environment or bath are 
an important tool. 

The study of master equations with time-dependent Hamiltonians is certainly not a new topic, going back at least as far as 
the pioneering work of Davies & Spohn lfT3ll . who derived an exact master equation for an adiabatic but infinitesimally weak 
system-bath interaction. Other, more recent approaches have been attempted, but in each case certain limitations apply. For 
example, Childs et al. used the Lindblad equation with a time-independent Hamiltonian at each time step as an approximation 
to the adiabatic evolution of a system with a time-dependent Hamiltonian lfT4l . Sarandy & Lidar derived a phenomenological 
adiabatic master equation, based on the idea that in the adiabatic limit the dynamical superoperator can be decomposed in 
terms of independently evolving Jordan blocks |15|. This approach is phenomenological in the sense that it does not allow 
one to derive the various parameters appearing in the final master equation from underlying system and bath Hamiltonians. 
Approaches based on non-Hermitian effective Hamiltonians, e.g., Refs. |fl6HT8ll , are necessarily also phenomenological. A 
rigorous phenomenological master equation derivation was given by Joye [ 19 1. Oreshkov & Calsamiglia connected open system 
adiabaticity to the theory of noiseless subsystems [20]. Thunstrom et al. derived a master equation from first principles in the 
physically reasonable joint limit of slow change and weak open system disturbances, but did not elucidate the relative time 
and energy scales involved in their approximations ETI . Various authors provided derivations for slow periodic Hamiltonians 
Il22lj24l . Such derivations are valuable and can be made rigorous, but the assumption of periodicity can be excessive, especially 
in the context of adiabatic quantum computation. Bounds on the validity of the adiabatic approximation for open systems, but 
without master equations, were presented in Ref. [25 1 (see also Ref. [10]). Various authors derived or studied adiabatic master 
equations limited to the case of a single qubit, where detailed physical considerations are possible I26 H281 . 

Our goal in this work is to derive master equations for adiabatic open system dynamics from first principles, while keeping 
track of all physical approximations, time- and energy-scales. In this manner we hope to fill a gap in the earlier literature on 
this topic, and to provide tools allowing for detailed comparisons with experiments satisfying the explicit assumptions behind 
our approximations. In particular, we derive several Markovian master equations, distinguished by different levels of adiabatic 
perturbation theory. When we add the rotating wave approximation (sometimes called the secular approximation) we arrive at 
master equations in Lindblad form [29], for which positivity of the state is guaranteed at all times [30|. Our formalism allows 
for the calculation of non-adiabatic corrections, which we also discuss. We apply our master equations to numerically study the 
evolution of a transverse field Ising chain coupled to a thermal bath, and discuss generic features of the evolution. 

Our basic starting point is the observation that the Markovian and adiabatic limits are fundamentally compatible, in the sense 
that a Markovian bath is "fast", while an adiabatically evolving system is "slow". As long as the corresponding timescales 
are appropriately matched, it is possible to derive an adiabatic master equation which is internally consistent. Somewhat more 
technically, we observe that in the interaction picture with respect to the unperturbed system and bath evolutions, where the 
bath is evolving "fast" while the system is evolving "slowly", and for sufficiently weak system-bath coupling, it is possible to 
consistently apply adiabatic perturbation theory to the time-dependent system operators. This insight allows us to derive our 
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adiabatic Markovian master equations. Our work is conceptually similar in its starting point to that of Amin et al. ||3"T1 (see 
also the Supplementary Information of Ref. Ell'), but is significantly more general. Some of our results are also conceptually 
similar to those found by Kamleitner & Shnirman [24 1 by de Vega et al. Il32"ll , and by Yung |33|, but are again more general. 
The master equations we derive are a natural generalization of standard time-independent Hamiltonian results 1 30 . 34 1, and our 
master equations reduce to the standard results after freezing the time dependence of the system Hamiltonian. 

The structure of this paper is the following. We set the stage in Section [II] where we define the system-bath model, perform 
the Born-Markov approximation, and introduce the bath correlation functions and spectral-density. We pay special attention to 
constraints imposed (via the KMS condition) by baths in thermal equilibrium, and single out the case of the Ohmic oscillator bath. 
We interrupt the formal development in Section [III] where we provide a summary of all the time- and energy-scales appearing 
in our various approximations, and the inequalities they must mutually satisfy. We then proceed to derive our adiabatic master 



equations in Section IV We proceed in several steps. First, in subsection IV A we find the adiabatic limit of the time-dependent 



system operators. Next, in subsection IV B we find two pairs of master equations in the adiabatic limit, one pair in the interaction 
picture, the other in the Schrodinger picture. The master equations within a given picture are distinguished by whether we apply 



the adiabatic approximation once or twice. Next, in subsection IV C we introduce the rotating wave approximation, which allows 
us to reduce the Schrodinger picture master equations into Lindblad form. We conclude the formal development by discussing 



non-adiabatic corrections to our master equations in subsection IV D We move on to a detailed numerical study of an example in 
Section[V] of a ferromagnetic chain in a transverse field, coupled to an Ohmic oscillator bath at finite temperature. We apply two 
of our Schrodinger picture master equations, with and without the rotating wave approximation, and show that without inclusion 
of the Lamb shift term they yield similar predictions. When the Lamb shift is included, however, substantial differences emerge. 
We discern four distinct phases in the evolution from the transverse field to the ferromagnetic Ising model, which we discuss and 



analyze. Concluding remarks are presented in Section VI The paper is supplemented with detailed appendixes where many of 
the technical details of the derivations and required background are provided, both for ease of flow of the general presentation 
and for completeness. 



II. PRELIMINARIES 



Model 



We consider a general system-bath Hamiltonian 



H{t) = H s (t) + H B + Hj 



(1) 



where Hs(t) is the time-dependent, adiabatic system Hamiltonian, Hb is the bath Hamiltonian, and Hi is the interaction 
Hamiltonian. Without loss of generality the interaction Hamiltonian can be written in the form: 



Hi = g^A a ®B a 



(2) 



where the operators A a and B a are Hermitian and dimensionless with unit operator norm, and g is the (weak) system-bath 
coupling. The joint system-bath density matrix p(t) satisfies the Schrodinger equation p = —i[H, p(t)], where we have assumed 
units such that H= 1. 

Let Us(t,t') = T + exp[— i f^dr Hs(t)] and UB{t,t') — exp(—i(t — t')Hs) denote the free system and bath unitary 
evolution operators (T + denotes time ordering), and define Uo(t,t') = U~s(t,t') ® UB(t,t'). Likewise, let U(t,t') — 
T + exp[ — i J dr H (r)] denote the joint Schrodinger picture system-bath unitary evolution operator. We transform to the 
interaction picture with respect to Hs(t) and Hb, by defining U(t, 0) = Uq (t, 0)U(t, 0), which, together with the interaction 
picture density matrix p(t) = Ul{t, 0)p(t)Uo(t, 0) satisfies 



^-U(t, 0) = -%Hi(t)U(t, 0) , U(0, 0) = 1 
at 



dt 



p(t) = -i[H I (t),p(t)}, p(0) = p(0) 



(3a) 
(3b) 



We restrict the use of tilde variables to refer to variables in the interaction picture. The time-dependent interaction picture 
Hamiltonian Hi(t) is related to its Schrodinger picture counterpart via: 



H T (t) = [rj(f,0)fl z [/o(t,0) - 9^2 Mt) ® B a{t) 



(4) 



3 



where we have defined the time-dependent system and bath operators: 

A a (t) = Ul(t,0)A a U s (t,0) 
B a (t) = U B (t,0)B a U B (t,0). 

We are interested in deriving a master equation for the system-only state, 

p s (t) = Tt B \p(t)] = Ul(t,0)Tr B (U B (t,0)p(t)U B (t,0))U s (t,0) 

= Ul(t,0)Tr B (U B {t,0)U B (t,0)p(t))U s (t,0) = ul(t,0)p s (t)U s (t,0) , 

where in the second line we used the fact that U B acts only on the bath. 



(5a) 
(5b) 

(6a) 
(6b) 



B. Born-Markov approximation 



Writing the formal solution of Eq. Ob} as 



p(t)=p(0)-i / drlH^rpir)), 



(7) 



and substituting this solution back into Eq. pb| ), we obtain, after tracing over the bath degrees of freedom, the equation of motion 
for the system density matrix ps{t) — Tr B p(t) 



-iTrr 



Hi(t),p(0) 



Tr, 



(8) 



We make the standard Born approximation assumption, that we can decompose the density matrix as p(t) — ps(t) <8> p B + x(i) 
where x(*)> which expresses correlations between the system and the bath, is small in an appropriate sense and can hence be 
neglected from now on [30 , 35, 36 1. Thus the equation of motion reduces to: 



(9) 



-p s (i) =g 2 V / dT(Af)(t-T)ps(t-T)A a (t)- A a (t)Ap(t-T)p s (t-T))B a p{t,t-T) + h.c. , 
at *,fi Jo 

where we defined the two-point correlation functions: 

B aP (t,t-T) = (B a (t)Bp(t-T)) = Tr[B a (t)Bp(t-T)p B ] ■ 



(10) 



In Eq. |9j, we have assumed for simplicity that (B a )o = Tr [B a (t)p B (Q)] — 0, so that the inhomogeneous term in Eq. ([8]) 
vanishes. Since we assumed that the bath state p B is stationary, the correlation function is homogenous in time: 



Bap(t, t - t) = {B a {t)B p {t - t)) = (B a (T)B (O)) = (B p {0)B a ( T ))* = B a0 {r, 0) 



(11) 



For notational simplicity we will denote S Q| g(r, 0) by B a p(r) when this does not lead to confusion. Let us denote the time scale 
over which the two-point correlations of the bath decay by t b , e.g., |B a ^(r)| ~ exp(— t/t b ). More precisely, we shall require 
throughout that 



POO 

/ dTT n \B afi (T)\~T n B + \ nG {0,1,2}. 



(12) 



As we show in Appendix [b] if t b <C 1/g, then we can apply the Markov approximation to each of the four summands in 
Eq. Q, i.e., replace ps(t — r) by ps(t): 



dT(...p s (t-T)...)B aP {r) 



dr{...p s {t)...)B atj {r) + 0{r B g z ), 



(13) 



where (...) refers to the products of A a and Ap operators in Eq. |9]), and where we have also assumed that t S> t b and the 
integrand vanishes sufficiently fast for t > Tg, so that the upper limit can be taken to infinity. Note that by Eq. ( p~2| > the integral 
on the RHS of Eq. ( fT3j ) is of 0(t b ), so that the relative magnitude of the two terms is 0[(gr B ) 2 ]. An explicit derivation of 
the upper bound on the error due to this approximation can be found in Appendix [B] The resulting Markovian equation cannot 
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resolve the dynamics over a time scale shorter than r B . 



C. Correlation functions, the KMS condition, and the spectral-density matrix 



In computing the terms appearing in Eq. (\3\ we are faced with integrals of the form J °° drAp(t ~~ T)ps(t)A a (t)B a p(T) and 
J °° d,TA a (t)Ap(t — T)ps(t)B a p{j). Our goal is to express these integrals in terms of the spectral-density matrix 



(14) 



r a/? H = / dre^B^r) 



the standard quantity in master equations. It is convenient to replace the one-sided Fourier transform in the spectral-density 
matrix by a complete Fourier transform. Thus we split it into a sum of Hermitian matrices, 



where we show in appendix [C]that 7 and S are given by 



dw' 
~2k 



(15) 

(16a) 
(16b) 



If we assume not only that the bath state is stationary, but that it is also in thermal equilibrium at inverse temperature /?, i.e., 
Pb = e~P HB J Z, then it follows that the correlation function satisfies the KMS (Kubo-Martin-Schwinger) condition l30l 



(B a (j)B b {0)) = (B b (0)B a (T + i/3)) 



(17) 



If in addition the correlation function is analytic in the strip between r = — i/3 and r = 0, then it follows that the Fourier 
transform of the bath correlation function satisfies the detailed balance condition 



(18) 



For a proof see Appendix [D] 

It is important to note that the KMS detailed balance condition imposes a restriction on the decay of the correlation function. 
To see this, note first that \B ab (-r) \ = ^Tr[p B U B (T)B a U B (T)B b ]\ = ^Tr[B a U B {T)B b U B (T)p B }\ = \{B a (0)B b ( T ))\ = 
\(B b (T)B a (0))\ = \B b a(r)\, where we used Eq. (JTTJ. Now assume that Eq. ( p~2] > would have to hold for all n. It would follow 



that 



d n 
duj n 



j=0 



r n e^B ab (r)dr\ u=0 



— OO 

n+1 



T n B ab (T)dr < / r n (\B ba (r)\ + \B ab (T)\) 

-00 Jo 



dT 



2r™ +1 VnG{0,l,...}, 



Thus all derivatives of 7 a t,(w) would have to be finite at u> = 0. 
On the other hand, it follows from the KMS condition ( p~8] > that 



d(—oj 



uj>0 



du> 



j>0 



so that in the limit as uj approaches zero from below or above, 

7 / a6 (0_)=/3 7 6a(0)-76a(0+)- 

This shows that in principle 7 a a (w) may be discontinuous at oj = 0. Indeed, the continuity condition 7^ o (0_) 
implies, from the KMS condition recast as Eq. ( |2T| ), that 

2^(0) = haa{0) • 



(19) 

(20) 
(21) 

7aa(0+) 

(22) 



This conclusion can be extended to the entire 7 matrix by diagonalizing it first. When Eq. p2| is not satisfied 7" a (0) diverges, 
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so that Eq. ( fT9| does not hold except for n £ {0, 1}. A simple example of this is j a b(u) — c > (constant) for oj > 0. 
Another example is a super-Ohmic spectral-density 7 a b(w) = w 2 /(l — e~ /3 ") for w > and 7 a b(w) = uj 2 /(e~ ,3u — 1) 
for a; < 0. Both examples violate Eq. ( fT9) ) for n > 2. Note, moreover, that when this happens, it follows from Eq. ( fT9] ) that 
Jo°° 1-2 (l^ab(~ T )l + l^af>( T )l) ^ T i s divergent, so that we must conclude that |i? a b (t) | > 1/t 3 for sufficiently large |t|, meaning 
that the correlation function has a power-law tail and, in particular, cannot be exponentially decaying. 

On the other hand, Eq. ( |22| > tells us that continuity and a lack of divergence at u> = require 7q O (0) to satisfy a condition 



which prohibits it from being arbitrary. This is indeed the case for the Ohmic oscillator bath discussed in subsection V A which 



satisfies Eq. ( |22| i with finite 7 QQ (0). For this case the bath correlation function is exponentially decaying assuming the oscillator 



bath has an infinite cutoff. However, as we show in subsection V B the Ohmic bath transitions from exponential decay to a 
power-law tail for any finite value of the cutoff, at some finite transition time r tr . In this case we find \B a p(r)\ ~ (t/tm) -2 , 
where tm is a time-scale associated with the onset of non-Markovian effects, and hence we have to relax Eq. ( fT2| ), and, replace 
it with 

11 dr T n \B aP (T)\ ~t£ +1 , n£{0, 1,...} (23a) 



/>oo poo 

/ dr\B afl (T)\~ / dT(r/T M )- 2 = T 2 M /r tI 

J Ttr J Tt.r 



(23b) 



III. TIMESCALES 

In this subsection we summarize, for convenience, the relations between the different timescales which shall arise in our 
derivation. The total evolution time is denoted tf and the minimum ground state energy gap of H$ is A, i.e., 

A= min ei(t)-e (t), (24) 
te[o,t/] 

where e {t) and E\ (t) are the ground and first excited state energies of Hg (t) . We shall arrive at master equations of the following 
general form: 

Ps(t) = [Ani(t) + £l(t) + CZT d (t)]Ps(t) ■ (25) 

■*ad 



Here C un [ = —i[Hs(t)+His-, •] is the unitary evolution superoperator including the Lamb shift correction, £jjf ss is the dissipative 

non-E 
diss 

Let 



superoperator in the fully adiabatic limit, and £j° s "" ad is the dissipative superoperator with leading order non-adiabatic corrections 



h= max \{e a (s)\d s H s (s)\e b (s))\ , (26) 

se[0,l];a,b 

where s = t/tf is the dimensionless time. To ensure that £ un ; generates adiabatic evolution to leading order we shall require the 
standard adiabatic condition 

A- f « 1 • (27) 

In order to derive Eq. ( pB} , the three superoperators are ordered in perturbation theory, in the sense that 

IIAndll » ll^dfjl » \\£Z*\\ , (28) 

where the norm could be chosen as the supoperator norm, i.e., the largest singular value (see Appendix [A|. We may then assume 
that 

IIAoniH » A » ||£*J . (29) 

Combining the 0(tb) due to the integral on the RHS of Eq. ( fT3] > (as already remarked there) with the g 2 prefactor from Eq. |9|, 
we have 

ll^dfssll ~ 9 2 tb ■ (30) 



6 



To ensure the first inequality in Eq. d28]l we thus require 



2 

^ « 1 ■ (31) 

The non-adiabatic correction is of order A ^ t , and when it appears in £J]° s "" ad it is multiplied by the same factor as £^ ss , i.e., we 
have 

ii£r; ad ii~s 2 r BA ^. 02) 

To ensure the second inequality in Eq. ( |28] i thus amounts to the adiabatic condition, Eq. ( p7| . All this is added to the condition 

gr B < 1 , (33) 

for the validity of the Markovian approximation, mentioned in the context of Eq. ( fT3j ), and justified rigorously in Appendix [B] 

There is one additional, independent time scale we have to concern ourselves with. This is the time scale associated with 
changes in the instantaneous energy eigenbasis relative to tb- If we require the change in the basis to be small on the time scale 
of the bath tb, we must have: 

-P^ « 1 . (34) 

We justify this claim in subsection IV A and AppendixlF] 

Note that the adiabatic condition, Eq. p7) , implies <C Atb, while Eq. p4| can be written as <C A -L. Putting this 
together thus yields 

^«mm(Ar B ,-l). (35) 
Atf Atb 

Our other inequalities [Eqs. ( f3T| , ( |33] >] can be summarized as 

gr B <min(l,A/p). (36) 



IV. DERIVATION OF ADIABATIC MASTER EQUATIONS 

A. Adiabatic limit of the time-dependent system operators 

Let us first work in the strict adiabatic limit. We will discuss non-adiabatic corrections in subsection IIVDI We denote the 
instantaneous eigenbasis of Hs(t) by {|e Q (i))}, with corresponding real eigenvalues (energies) {s a (t)}, i.e., Hs(t) |er a (t)) = 
£a(t) \e a (t)), and Bohr frequencies w& Q (i) = £b(t) — s a {t). As shown in Appendix E 1 we can then write the system time 
evolution operator as: 



Us (t,t') = Uf(t,t') + 0^-^-^ , 



(37a) 



Uf(t,t') = E k.(*)><e-(0|e- Wt, ° . (37b) 

a 

where Ug d (t,t') is the "ideal" adiabatic evolution operator. It represents a transformation of instantaneous eigenstate \e a (t')) 
into the later eigenstate \e a (t)), along with a phase 

MM') = / dT[e o (r)-0„(T)] , (38) 



E2 



where (j) a {t) — i(e a (t)\e a {t)) is the Berry connection. The correction term O (^ AH f ) ^ s derived in Appendix '. 

To achieve our goal of arriving at a master equation expressed in terms of the spectral-density matrix, our basic strategy is to 
replace the system operator Ap(t — r) = Ul(t — r, 0)ApUs(t — r, 0) with an appropriate adiabatic approximation, which will 
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allow us to take this operator outside of the integral. To see how, note first that 

U s (t - r, 0) = U s (t - r, t)U s (t, 0) = Ul{t,t - T )U 8 {t, 0) . (39) 

We now make two approximations: first, as per Eq. ( |37a[ ) we replace Us(t, 0) by Ug d (t, 0); second, we replace Ug(t, t — r) by 
e lTHs w, an approximation justified by the appearance of the short-lived bath correlation function B a0 (r) inside the integrals 
we are concerned with. Thus, we write 

U 8 (t-T,0) = e lTHs ^Uf(t,0) + Q(t,T) , (40) 
and find the bound on the error due to dropping Q(t, r) in Appendix|F] Let 

H ba (t,t') = Vb(t,t') - Ha(t,t') , (41a) 
U ab (t) = \e a (t))(e b (t)\. (41b) 

Neglecting the operator-valued correction term Q(t, r) entirely, we have, upon substituting Eq. ([37b) and using e lrHs w \e a (t)) = 
e ire »W|e a (t)),that 

dTA (t-T)p s (t)A a (t)B a0 (T) (42a) 
dTUf(t,0)e- lTHs ^Ape iTHs ^Uf(t 7 0)p s (t)A a (t)B aP (T) (42b) 
= / dTY^e-^^eamiEaWe-""^^^"^ (42c) 

J ° ab 

= ^e- Vb " (t ' 0) |ea(0))(e a Wl^kfcW)(£6(0)|ps(t)A Q (i) / dre iT ^- e ^B a0 (r) , (42d) 







ab 



where the approximation in ( |42b| > is shown in Appendix F to be O [tb min{ 1 , + ^f^ }] ■ The first term, -^i^ , is the smallness 

t 2 h 

parameter of the adiabatic approximation, which we have already assumed to be small. The second term, -f— [mentioned in 
Eq. (134)], is new, and its smallness is associated with changes in the instantaneous energy eigenbasis relative to tb- We are 
interested in the regime where both terms are small. 
Thus, to the same level of approximation 

/>oo 

/ drA p {t ~ T)p s (t)A a (t)B a p(T) « V e'^^A^n^ps^A^r^iubait)) , (43) 

J ° ab 

where 

A aab (t) = (e a (t)\A a \e b (t)) = A* aba (t) , (44) 
and T a p(u>b a (t)) is the spectral-density matrix defined in Eq. ( fl4) . Similarly, we have for the other term 



/>oo 

/ drA a (t)Ap(t - T)~ Ps (t)B alj {T) » V e-^^A^A^Uab^ps^T^bait)). 



(45) 



B. Master equations in the adiabatic limit 

We are now ready to put everything together. Starting from the Born-Markov master equation constructed from Eqs. |9) 
and ( p"3) , and using the approximations ( |43) and ( |43) , we arrive at the following one-sided adiabatic interaction picture master 
equation: 

j t ~p s {t) = g 2 J2 e-^ (t ' 0) J2 r a pMt))A Pab (t) [U ab (0)p s (t),A a (t)} + h.c. . (46) 

ab a.j3 

Since we used an adiabatic approximation for Ap(t — r), it makes sense to do the same for A a (t), i.e., to replace the latter 
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with (t, 0) A a Ug d (t, 0). If this is done, we obtain the double-sided adiabatic interaction picture master equation 

j t Ps{t) =g 2 Y, e-^-f*' '-^*' )] Y, T <*pMt))A acd (t)Ap ab (t) [U ab (0)p s (t), n od (0)] + h.c. . (47) 

abed a/3 

It is convenient to transform back into the Schrodinger picture. Using ps(t) = Ug(t, 0)ps{t)Us{t, 0) [Eq. |6|] implies that 

dlPsif) = Ug(t, 0) (i[Hs(t), ps(t)] + jliPs(t)) Us(t, 0). Hence, using Eq. ( |46| ), we find the one-sided adiabatic Schrodinger 
picture master equation 

j t p s {t) = -i[H s (t),ps{t)} +9 2 YY, r ^KW) [Jw(*)flsr(t)> M + h.c. , (48) 

ab a{3 

where 

CpW = e-^WApatifiUsit, 0)n ab (0)ul(t, 0). (49) 

This form of the master equation has not appeared in previous studies of adiabatic master equations. If we again use the 
adiabatic approximation for U s (t, 0), i.e., replace U s (t, 0)IL ab (0)Ul(t, 0) by Uf*(t, 0)Il ab (0)Uf(t, 0), we obtain the double- 
sided adiabatic Schrodinger picture master equation 

j t Ps{t) = -% [H s (t),p s (t)} +9 2 J2J2 r «tf(w&a(*)) [L abt0 (t)ps(t),A a ] + h.c. , (50) 

a/3 ab 

where 

L ab<a (t) = A aab {t) \e a {t))(e b (t)\ = L\ a a {t) . (51) 

Comparing to Eq. ( |25j ), the first term in Eqs. ( |4"8j ) and ( [50] ) is C ua i(t), while the second is £^ ss (t). 

The master equations we have found so far are not in Lindblad form, and hence do not guarantee the preservation of positivity 
of ps- We thus introduce an additional approximation, which will transform the master equations into completely positive form. 



C. Master equation in the adiabatic limit with rotating wave approximation: Lindblad form 

In order to arrive at a master equation in Lindblad form, we can perform a secular, or rotating wave approximation (RWA). 
To do so, let us revisit the t oo limit taken in the Markov approximation in Eq. ([13]). Supposing we do not take this limit 
quite yet, we can follow the same arguments leading to Eq. ( |42c| i, which we rewrite, along with the adiabatic approximation 
A a {t) « Uf\t, 0)A a Uf(t, 0). This yields 

t 

drA p {t - T)p s {t)A a {t)B a0 {T) a (52a) 
f rf^ 5^ e-*^- C £ '°> ° | &a (0)> <e„ (t) ] j ff6 (t)> <e b (0) 1 p s (t) | &c (0) ><e c (t) | ^« [ e rf (t) ><e d (Q) 1 e^™^ S a/3 (^) . 

•'° abed 

(52b) 

We note that p dc (t, 0) + p. ba (t, 0) = /„ dr [w dc (r) + u> ba (r) - (4> d (r) - 4> c {t)) + (0 h (r) - (j) a (r))]. One can now make the 
argument that when the f oo limit is taken, terms for which the integrand vanishes will dominate, thus enforcing the "energy 
conservation" condition tu ba = —LUd c - This is a similar rotating wave approximation as made in the standard time-independent 
treatment, although here, the approximation of phase cancellation is made along the entire time evolution of the instantaneous 
energy eigenstates. Clearly, in light of the appearance of other terms involving t in Eq. ( |52b| i, this argument is far from rigorous. 
Nevertheless, we proceed from Eq. (|52b[) to write, in the t — > oo limit, 



drA fi (t - T)p s (t)A a (t)B a0 {T) » Y ^, u W^^(<)nu,(0)ps(t)n u (0)r^(w) , (53) 



9 



where we have defined a new index u> such that: 

A a , u (t)= Yl (e a (t)\A a \e b (t)) , 1^(0)= £ |e Q (0))(£ b (0)| . (54) 

£&(*)— £a(t)=u: et,(t) — e a {t)=u: 

Note that the set of {w}'s involved in the sum i s evolving in time since it corresponds to differences of the instantaneous 
energy eigenvalues, but we suppress the time dependence for notational brevity. We show in Appendix |G|how, by performing a 
transformation back to the Schrodinger picture, along with a double-sided adiabatic approximation, we arrive from Eq. ( |53] l at 
the Schrodinger picture adiabatic master equation in Lindblad form: 

p s (t) = -i [H s (t) + H hS (t) , ps (t)] + £ £ lap (w) 

a/3 w 

where the Hermitian Lamb shift term is 
and we have defined 

L u , a (t)= LabAt)- (57) 

e b (t)-e a (t)=ui 

Since the bath correlations are of positive type, it follows from Bochner's theorem that the matrix 7 — the Fourier transform 
of the bath correlation functions — is also positive QUI . Therefore, this Lindblad form for our master equation guarantees the 
positivity of the density matrix. 

We emphasize that Eqs. ( |48j ), ( |50| l, and ( |55] l all generalize both the standard Redfield and Lindblad time-independent Hamil- 
tonian results to the case of a time-dependent Hamiltonian in the adiabatic limitQ The time-independent result can easily be 
recovered by simply freezing the time dependence of the Hamiltonian, energy eigenvalues and eigenvectors. To see this explic- 
itly, let us restrict ourselves to the Lindblad case. 

The energy eigenvalues and eigenvectors are time independent in this case, so we can replace our time dependent Lindblad 
operators with time independent ones, and the ui index no longer varies with time: 

L u , a (t) -> L u , a = \ £ a)(£a\A a \e b )(e b \ ■ (58) 

The resulting equation becomes: 

,p,ps(t)} , (59) 

a/3 ui ^ ' 

which is the standard form for the time-independent Lindblad master equation. This should make evident the physical meaning 
of our derivation. We have systematically generalized the time-independent result such that the Lindblad operators now rotate 
with the (adiabatically) changing energy eigenstates, which makes them time-dependent. This is a non-trivial difference, as it 
encodes an important physical effect: the dissipation/decoherence of the system occurs in the instantaneous energy eigenbasis. 

Equations (|50) and ( |5"5] ) are the two master equations we use for numerical simulations presented later in this paper. We note 
that Eq. ( |55| l appears similar to the Markovian adiabatic master equation found in Ref. E4ll . but is more general and did not 
require the assumption of periodic driving. 



L Ut0 (t)ps(t)Ll ta (t) - \ {Li, a (t)L Ui p(t),p s (t)} 



(55) 



D. Non-Adiabatic Corrections to the Master Equations 



So far, we assumed the adiabatic limit of evolution of the system [see Eq. ( |37a| i]. The adiabatic perturbation theory we review 
in Appendix |E 1| allows us to compute systematic non-adiabatic corrections to the master equations we have derived. This 



perturbation theory is essentially an expansion in powers of and we rederive in Appendix E 1 the well known result [5] 



1 See for example Eq. (3.143) in Ref. (30) and Eq. (8.1.33) in Ref. (34). 



10 



that to first order 

U s (t,t') = Ug(t,t') [I + Vi(t,0] > (60) 

where 



^ 1 (t,<') = -Ei^( < '))^^)i/, tdre " ,;AIbQ(T,0 ( e -( r )i^( r )) 



(61) 



Thus, to de rive th e lo west o rder non-adiabatic corrections to our master equations is a matter of repeating our calculations of 



subsections 



IV A 



and 



IV C with Ug d (t,t') replaced everywhere by the first order term Ug d (t,t')Vi(t, t'), and adding the result 
to the zeroth order master equations we have already derived. Rather than actually computing these corrections, let us estimate 
when they are important. 



The condition under which the zeroth order adiabatic approximation is accurate is Eq. ( |27| i, which is now replaced by 

r < A 2 , (62) 



i.e., with the -C replaced by a mere <. However, we would still like to perform the approximation of Eq. ( |40| ), in the sense that 
U s (t - r,0) « e lTiIs ^Uf{t, t') [I + Vi(t, t')]. This still requires 

h 1 

r < -a , (63) 

which is what allows the use of e iTff s(*) m thj s approximation (as shown in Appendix]?]). Taken together, these two conditions 
are weaker than Eq. ( (35] l, which we can rewrite as jj <C min(A 2 , ^-). 

Recall that to ensure the validity of the Markov approximation and ||£ un i|| ^ ll^dissll' we a ^ so t0 demand the inequalities 
in Eqs. ( [3~Tj ), ( [33] >; these can be now summarized as max (g, <C to be compared with Eq. ( |63] l. 



V. AN ILLUSTRATIVE EXAMPLE: TRANSVERSE FIELD ISING CHAIN COUPLED TO A BOSON BATH 

A. The Model 

We proceed to use our master equations to study the evolution of the Ising Hamiltonian with transverse field 

H s (t) = A(t)H§ + B(t)H§ , (64a) 

N 

#* = ^>f, (64b) 



i=l 

N N 

where the functions A(t) and B(t) are shown in Fig.[T] and were chosen for concreteness to describe the interpolation between 
the transverse field and Ising term in the D-Wave Rainier chip [2|. This is a system which begins with the transverse magnetic 
field Hg turned on while the Ising term Hg is turned off, and then slowly switches between the two. 
We couple this spin-system to a bath of harmonic oscillators, with bath and interaction Hamiltonian 

oo N 

fc=l f=l k 

where bt and are, respectively, raising and lowering operators for the fcth oscillator with natural frequency uif., and g J k is the 
corresponding coupling strength to spin j. This is the standard pure dephasing spin-boson model [37 1, except that our system is 
time-dependent. The resulting form for our operator L [Eq. |5T)] is 

(t)){s a (t)\a!\e b (t)){e b (t)\ = A tab (t)\e a (t))(e b (t)\ . (66) 
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Recall that our analysis assumed that the bath is in thermal equilibrium at inverse temperature /? = l/(fcgT), and hence is 
described by a thermal Gibbs state ps = exp (— {3Hb) jZ. We show in Appendix |H| that this yields 

T«M = ^St*I*I ( e <<"> + «-" M e(-«)) (67a) 

s «<*«» - f <^ot" ( p (dr^) + e "" p (dk)) • <67b) 

where only one of the Heaviside functions is non-zero at to = 0. To complete the model specification, we assume an Ohmic bath 
spectral function 

J(cj) = rjuje'^/^ , (68) 

where uj c is a high-frequency cut-off and rj is a positive constant with dimensions of time squared. 

It is often stated that the terms associated with the Lamb shift H^s [Eq. (|56]l], i.e., Eq. (|67b|i, can be neglected, since the 



relative order of S and Hs is (; 2 tb/A, and indeed we have assumed g 2 rs/A <C 1 [Eq. (31 1], However, this argument is 
misleading for two reasons. First, S can be divergent, as is easy to see in the limit u c = oo for the Ohmic model |68|, where 
for uj 3> maxba ojj, a (t), the integrand tends to a constant. Second, S should be compared to 7, as both are of the same order 
g 2 Ts, and both result in changes to the system relative to its unperturbed state. Indeed, in the interaction picture with respect to 
Hs + -Hls [recall Eq. (|55j)], the overall transition rates between states with quantum numbers a and b will depend on the dressed 
(i.e., shifted) energy gaps lj^ + tjff. The importance of this Lamb shift effect was also stressed by de Vega et al. (32). We 



analyze the Lamb shift effect in subsection V C 



Finally, we note that although the harmonic oscillators bath with linear coupling to the system provides a 7 matrix that 
satisfies the KMS condition, it is important to note that this model has infrared singularities that destroy the ground state of the 
total system 1 38 1. The KMS condition assumes a stable ground state and stable thermal states, which our underlying spin-boson 
model violates. However, for the purposes of our work, a 7 matrix that satisfies the KMS condition will suffice without too much 
concern about how it is derived. 



B. Correlation function analysis 

In light of the subtleties alluded to in subsection |II C| associated with satisfying the KMS condition, we analyze the different 
timescales determining the behavior of the Ohmic correlation function in this subsection. Removing the oj dependence from the 
g's in Eq. ( |67a| i and substituting J(oj) from Eq. ( |68] l, it is possible to compute the bath correlation function analytically for the 
resulting 

- ™-e-?» ffV ' (69) 

by inverse Fourier transform of Eq. ( |16a[ ). The result is 

B ab (r) = ( ^ (^- + %)+ (l + J- ~ %)) , (70) 
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where tjA" 1 ^ is the mth Polygamma function (see Appendix [i] for the derivation). We first assume 

/3lo c » 1 , (71) 

and consider an expansion in large f3uj c : 

- {-^' (t) + | *T-5<^3 ■ 

followed by an expansion in large t/(3: 



B ab (r) = {-^e-^» + + O (e- 2 ^,r-3)) . (73) 

This expansion reveals the two independent time-scales that are relevant for us. First, there is the time scale tb associated with 
the exponential decay (corresponding to the true Markovian bath), given by: 

t b — > t b = — , (74) 

Z7T 

then the time scale associated with non-Markovian corrections (the power law tail): 

tm = J—. (75) 
V w c 

For sufficiently large uj c , these two time scales capture the two behaviors found in S (,(t), as illustrated in Fig.[2ja). 

The transition between the exponential decay and the power law tail occurs at a time r tr given by 4n 2 e~ Ttr / TB = (^) 2 , or 
equivalently by 

1® = ^., »-*r^. CO 

This transcendental equation has a formal solution in terms of the Lambert-VK function (39), i.e., the inverse function of f(W) = 
We w , as can be seen by changing variables to y = —9/n, and rewriting 9 n e~ e = a = 2/((3u> c ) as ye v = —a 1 /™/?!, whose 
solution is 9 = —ny = — nW(— ^a 1 /™). However, for our purposes the following observations will suffice. We seek a 
Markovian-like solution, where Tt r is large compared to the thermal timescale set by j3, i.e., we are interested in the regime 
where 9 ^> 1. In this case we can neglect 9 2 compared to e e , and approximate the solution to Eq. ( fTo} by 9 ~ ln(/?w c /2). Thus 

r tr ~ j81n(^o; c ) . (77) 

This agrees with the first term of the asymptotic expansion W(x) = ln(x) — lnln(a;) + 1 " I |"^' > + •••> which is accurate for 
x>3j39l. 

When f3w c ^> 1 is not strictly satisfied, the exponential regime is less pronounced, and t' b is corrected by powers of u> c . By 
dimensional analysis, the corrections must be of the form: 



where c n are constants of order one that must be fitted [see Fig.^b)]. 

The implications of this cutoff -induced transition for our perturbation theory inequalities are explored in Appendix [F] where 
we show that a sufficient condition for the theory to hold is 

' < m in{2r B ,^ + ^}, (79) 



w c ln(/3w c ) 1 A 2 tf tf 

which can be interpreted as saying that the cutoff should be the largest energy scale. Equation dWk joins the list of inequalities 



given in Section III as an additional special condition that applies in the Ohmic case, along with f3uj c ^> 1. 
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ln(S(r)) ln(S(f)) 




(a) u) c = 32ttGHz (b) cu c = SttGHz 



FIG. 2. (a) An example of the bath correlation function B a b(r) for /3 = 2.6 _1 ns and iv c — 327rGHz. The solid black curve is the function 
In \ B a t{r)\. The blue dashed curve is the function ln47r 2 exp(— 2tyt//3) and the red dot-dashed curve is the function In 2/?r -2 /uj c . (b) Same 
as in (a), but lo c — 87rGHz, and the blue dashed curve is the function In do exp(— (27t//3 + ciu; c + C2<^ 2 /3)t) with ci ~ —0.039, ci ~ —0.004, 
do « 33.13. 



C. Numerical Results 

For concreteness, we take — g, ui c = 87rGHz, and T = 20mK w 2.6GHz (in units such that H = 1; this is the operating 
temperature of the D-Wave Rainier chip [2Q, corresponding to tb = 0.06ns for the Ohmic model with infinite cutoff, Eq. ( f74] >, 
and tb ~ 0.07ns for oj c = 8tt using Eq. ( |78j ). For this value of u> c the transition between the exponential and power law regimes 
is still sharp (see Fig[2|b)) and occurs at approximately T tr = 0.33ns. For these values, we satisfy at least one of the cases from 
Eq. ((79]), including numerical prefactors: - r^fg^ s < 2tb- 

We focus on the N = 8 site ferromagnetic chain with parameters: 

h = - , h i>0 = , J i>i+1 = -1 , i = 0, . . . , 7 , (80) 

where we pin the first spin in order to break the degeneracy in the ground state of the classical Ising Hamiltonian. The system is 
initialized in the Gibbs state: 

e -PH s {0) 

Ps(t = 0) = . (81) 

To help the numerics, we truncate our system to the lowest n = 18 levels (out of 256), rotating the density matrix into the 
instantaneous energy eigenbasis at each time step. The error associated with this is small as long as our evolution does not 
cause scattering into higher n states, as we have checked. The forward propagation algorithm used is an implicit second order 
Runge-Kutte method called TR-BDF2 with an adaptive time step I40[ |4D . 

Figure [3]presents our results for the evolution of the system described in the previous subsection. We computed the overlap 
between the instantaneous ground state of Hs(t) and the instantaneous density matrix predicted by our two master equations ([50]) 
(non-RWA) and ( |55] ) (RWA). Although our two master equations predict different numerical values for this overlap, the qual- 
itative features of the evolution are the same. We observe a generic feature of four distinct regions of the evolution: a gapped 
phase (labeled "1" in figure [3}, an excitation phase (labeled "2"), a relaxation phase (labeled "3"), and finally a frozen phase 
(labeled "4"). We elaborate on these regions in the following subsection. Furthermore, we observe that the larger tf (therefore 
more adiabatic) evolution shows a smaller difference between the two master equations for the final fidelity. The results in Fig. [3] 
illustrate the importance of the Lamb shift. We find that while its effect is small in the RWA, its effect is relatively large in the 
non-RWA. 

In order to study the importance of the relaxation phase, we can study the behavior of the fidelity at t = tf as we change 
the coupling strength g. We observe (see Fig. |4| that there is a rapid drop in fidelity from the closed system result as soon as 
the coupling to the thermal bath is turned on, but there is a subsequent steady increase in the fidelity as the coupling strength 
is further increased. This increase in fidelity is a direct consequence of the higher importance of the relaxation phase (made 
possible by the increasing coupling strength) in restoring the probability of being in the ground state. However, we also observe 
that there is a very pronounced difference between the behavior of the results from the two master equations as the coupling 
is increased. In the case of the Lindblad equation, the fidelity saturates, whereas for the NRWA equation, we see an increase 
in fidelity and a subsequent violation of positivity. These results bring to light the relative advantages and disadvantages of 
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(e (t)\p(t)\s {t)) 
1.0 



0.2 0.4 0.6 0.8 

(a) Including Lamb shift 



t 

t.O T f 




0.2 0.4 0.6 0.8 

(b) Excluding Lamb shift 



t 

i.o T f 



FIG. 3. Fidelity, or overlap of the system density matrix with the instantaneous ground state along the time evolution for tf = 10/^s and 
rjg 2 /(h 2 ) — 1.2 x 10~ 4 /(2-7r). The solid blue curve was calculated using Eq. \50\ (no RWA), while the dashed red curve was calculated 
using Eq. {55} (Lindblad form, after the RWA). Four phases are indicated: thermal (1), excitation (2), relaxation (3), and frozen (4). Panel (a) 
includes the Lamb shift terms, while (b) excludes them. 



both master equations. For the Lindblad equation, positivity of the density matrix is guaranteed, but it clearly is not capturing 
the physics associated with the increasing importance of the thermal relaxation that the non-RWA equation captures. However, 
the non-RWA equation fails to preserve positivity of the density matrix so it is unable to reliably describe the system at higher 
coupling strength. Others have also noted that the RWA and NRWA can lead to physically different conclusions, e.g., in the 
context of Berry phases in cavity QED 1 421 . 

{e (tf)\p(tf)\Eo(tf)) (so(tf)\p(tf)\E (tf)) 




0.1 0.2 0.3 0.4 0.5 n i O'.OOOO 0.0001 0.0002 0.0003 0.0004 0.0005^2 



(a) RWA (b) non-RWA 



FIG. 4. Dependence of the fidelity at tf = 10/is on the coupling strength to the thermal bath is varied using our two master equations. The 
insets are closeups of the behavior near zero coupling. 



D. The Four Different Phases 

1. Phase 1 - the gapped phase 

For times sufficiently close to the initial time, the ground state of Hs(t) is the ground state of Hl$ , i.e., the state |0) = 
<3jLi \-)j, where |±) . = ± |t),-)/\/2 with energy e (0) = -NA(0), and where , are the +1, -1 eigenstates of 
a? (computational basis states of the jth spin or qubit). The lowest lying energy states are then the iV-fold degenerate states with 

a single flip of one of the spins in the x direction, i.e., \i) = |— ),• |+}j ®f=i+\ \~)j> w ith energy £i(0) = — (N — 2)^4(0). 
Therefore the gap between the ground state and the first excited states is: 

t<if 

A{t) =ei(t)-e (*) « [-{N- 2)A{t)]-[-NA(t)]=2A{t) , (82) 
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which is at least twice as large as our fc^T ~ 2.6GHz, almost until A(t c ) — B(t c ) « 5GHz at t c w 0.35t/ (see Pig. |TJ>. Noting 
that erf = | : F) i , we can write the following relations in terms of the ground and first excited states: 



uj i0 = A = -cj i , cr| |«) = |0) 



(83) 



Noting that erf \j) is a doubly excited state Vj 7^ i € {1, . . . , TV}, we truncate the problem to the ground and first excited states 
only, so that there are just three types of values of 7(0;) we need to consider: 7(0), 7(w<o)j li^Oi)- Recalling the KMS condition, 
Eq. ( p"8| ), we have 



7(^oi) 



7(^0) , 



(84) 



which shows that upward transitions are exponentially suppressed relative to downward transitions, by a factor ranging between 



-2A(0)/(k B T) 



-67.4/2.6 _ 7 x 10 -12 and e -2A(0.3i/)/2.6 



-3.8 



0.02. This explains why for early times (Phase 1) the 



system hardly deviates from the ground state, which in turn is very close to the thermal state dSTb. 




0.002 0.004 0.006 0.008 0.010 t f 

(a) Ground State 




t 

0.002 0.004 0.006 0.008 O.OlcT t f 

(b) Excited states 



FIG. 5. Early evolution of diagonal elements of the system density matrix in the thermal regime. Red dots are from using the Lindblad equation 
J55), while the solid blue curve is from the solution of Eq. l|85[>. 
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(b) Suppression factor 



1.0 tf 



FIG. 6. The behavior of the gap and the exponential suppression factor along the time evolution. The dashed line is ksT w 2.6GHz. 



To make this argument more precise, denoting poo = (0| p| 0) and pa = (i\p\i), we can write the effective (truncated to the 
ground and first excited states) Lindblad equation ( |5"5| l as the simplified rate equations 

Pit ~ 7w( w io) {-Pu + e~ PLUz0 Poo) , (85a) 
poo w 1 - N Pli , (85b) 
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where we have assumed that the system is initially in the thermal state (JHTJ and the gap is large (relative to fc^T). A derivation 
of Eq. ( |85| l can be found in Appendix |7] We compare our simulation results with the results from the above equations in Fig. [5] 
and find very good agreement for early times. 



2. Phase 2 — the excitation phase 

When A(t) becomes small enough such that /3A(t) ~ 0(1), then the KMS condition no longer suppresses excitations from 
the ground state to higher excited states (see Fig. [6}. If we interpret the master equation as a set of rate equations for the matrix 
elements of p, we can identify the rate of scattering into the ith state from the jth state as being the term in the pa equation with 
coefficient pjj. Therefore, we find that the rate of scattering from the ground state to excited states is given by 

Excitation rate oc 7(A(i))e~^ A( *Voo , (86) 

whereas the relaxation rate is given by 

Relaxation rate oc ^(A(t))pu . (87) 

Therefore, as we emerge from the gapped phase, we have p 0() ^> pa, so scattering into excited states dominates over relaxation 
into the ground state. This explains the loss of fidelity in Phase 2. 



3. Phase 3 - the relaxation phase 

As the gap begins to grow again and the suppression factor shrinks (see Fig. [6]), the KMS condition begins to suppress 
scattering into higher excited states while allowing relaxation to occur. In our model this causes a resurgence in the overlap with 
the ground state. Therefore, in this phase, the presence of the thermal bath can actually help to increase the fidelity, as was also 
observed in Ref. [|26l |32l |43l 

The excitation and relaxation phases reveal the two competing processes for a successful adiabatic computation. If we spend 
too long in the excitation phase (or if the gap shrinks too fast relative to tf and the evolution is not adiabatic), the system loses 
almost all fidelity with the ground state, and the system would have to spend a very long time in the relaxation phase to recover 
some of that fidelity. 

However, we stress that fidelity recovery will not be observed if the population becomes distributed over a large number of 
excited states in the excitation phase. This would happen, e.g., if when the gap closes there is an exponential number of states 
close to the ground state, such as in the quantum Ising chain with alternating sector interaction defects l44l . To see this more 
explicitly in the context of our analysis, note from Eq. ( |85b| i that if there is an exponential number N of pa coupled to poo (i-e., 
N is an exponentially large fraction of the dimension of the system Hilbert space), then poo decreases exponentially. By Eq. ([86]) 
this means that all pa become exponentially small, but not zero. In phase 3, by Eq. ((87), the relaxation is suppressed as long as 
the gap is not very large, because the relaxation is proportional to p ii7 which is exponentially small. This analysis presumes that 
the system-bath Hamiltonian has non-negligible coupling between the ground and excited state, i.e., that |(0|<r||i)| > in our 
model [see Eq. ( |J6a| l]. This suggests another mechanism that can suppress relaxation: the ground state in phase 1 might have a 
large Hamming distance from the ground state in phase 3. Relaxation is then suppressed simply because the coupling is small. 

We might expect that there exists an optimum tf for which the fidelity is maximized by the end of the relaxation phase. 
However, for the simple case of the spin-chain we considered, we did not observe such an optimum t f. The fidelity continues to 
grow (albeit slowly) for sufficiently high tf. This is illustrated in Fig. [7] 



4. Phase 4 - the frozen phase 

As the gap continues to grow, the relaxation phase ends [notice that the tail in Fig. |6|b) is longer than the actual re- 
laxation phase] and the system's dynamics are frozen in the ground state. This is because H$ becomes almost entirely 
diagonal in the a z basis, and so the off-diagonal components of the L^ i operators vanish (or become very small), i.e., 
Aiab(t) = ( e a(*)| <7 i \ £ b(t)) S a t [Eq. ( |66"] l] leaving only the diagonal ones. At this point, while off-diagonal elements may 
continue to decay, the system ground state is no longer affected by the bath; this is the onset of the frozen phase. 



17 



<e (Olp(0|eo(0) 
1.0 
0.8 
0.6 
0.4 
0.2 

0.2 0.4 0.6 0.8 1.0 t f 

FIG. 7. Time evolution of the system density matrix using the RWA equation with r/g 2 / (ft 2 ) = 0.4/ (2tt) for different tf's. Solid blue curve 
is for tf — 10/^s, dashed red curve is for tf = 100/is, and dot-dashed green curve is for tf = 1ms. 



E. Thermal equilibration 

An interesting question is whether the system reaches thermal equilibrium throughout the evolution. To answer this we 
computed the trace-norm distance [defined in Eq. < | A 1 1 >] between the instantaneous system density matrix and the instantaneous 
Gibbs state pabbs(t) = exp[— /3Hs(t)]/2>, where Z = Tr [exp (—j3Hs(t))} is the partition function, for the Ising chain discussed 
above. The result is shown in Fig. [8] The distance is zero at t = since the system is initialized in the Gibbs state, and then 
begins to grow slowly as the system transitions from the gapped phase to the excitation phase. Though not generic, the distance 
decreases as the gap shrinks while the excitation phase becomes the relaxation phase, and the system returns to a near Gibbs 
state where the gap is minimum (at t/tt ~ 0.4). As the gap opens up again in the transition from the relaxation phase to the 
frozen phase, the distance begins to grow and continues to grow throughout the frozen phase. As such, the system is quite 
far from the Gibbs state at the final time. This feature is significant for adiabatic quantum computation, since it shows the 
potential for preparation of states which are biased away from thermal equilibrium towards preferential occupation of the ground 
state. However, we note that on sufficiently long timescales one should expect (from general thermodynamic arguments) terms 

llpW-PGibbsWIli 

o.2o ; 

0.15; / 
o.io; / 

0.05 



FIG. 8. Trace-norm distance between the evolving system density matrix (using the Lindblad equation) and the Gibbs state for an annealing 
timeof^ ■ = lOfis and r;g 2 / (h 2 ) = 0.4/(2tt). 

proportional to a x and a v in the system-bath interaction, which we neglected in writing the interaction Hamiltonian Hj in 
Eq. ( |65| ), to become important, and to disrupt the frozen phase, allowing the system to fully equilibrate into the Gibbs state. 



VI. CONCLUSIONS 

Using a bottom-up, first principles approach, we have developed a number of Markovian master equations to describe the 
adiabatic evolution of a system with a time-dependent Hamiltonian, coupled to a bath in thermal equilibrium. Our master equa- 
tions systematically incorporate both time-dependent perturbation theory in the (weak) system-bath coupling g, and adiabatic 
perturbation theory in the inverse of the total evolution time tt. Since we have kept track of the various time- and energy-scales 
involved in our approximations, higher order corrections (starting at third order in g and second order in 1/tj) can be incorpo- 




0.2 
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rated if desired, a problem we leave for a future publication. Using two of our master equations, we studied generic features of 
the adiabatic evolution of a spin chain in the presence of a transverse magnetic field, and coupled to a bosonic heat bath. We 
identified four phases in this evolution, including a phase where thermal relaxation aids the fidelity of the adiabatic evolution. 
We hope that this work will prove useful in guiding ongoing experiments on adiabatic quantum information processing, and will 
serve to inspire the development of increasingly more accurate adiabatic master equations, going beyond the Markovian limit. 
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Appendix A: Norms and inequalities 

We provide a brief summary of norms and inequalities between them, as pertinent to our work. For more details see, e.g., 
Refs. 145H471 . Let \A\ = V A^A. Unitarily invariant norms are norms that satisfy, for all unitary U, V, and for any operator A: 
\\UAV \\ u i — ||-i4||ui- Examples of unitarily invariant norms are the trace norm 

\\A\\ 1 ^Tt\A\=J2^(A), (Al) 

i 

where Si(A) are the singular values (eigenvalues of |j4|), and the supoperator norm, which is the largest eigenvalue of |j4|: 

114100= sup J^\aTa\4>) =ma XSi (A) , (A2) 

Therefore \\A \ip) \\ < ||4loo for all normalized states \tJj), and ||-A||oo < WMi- 

All unitarily invariant norms satisfy submultiplicativity: 

||AB|| ui < IHIuiH-BlU ■ (A3) 

The norms of interest to us are also multiplicative over tensor products and obey an ordering: 

WA^BIU = ||4li||#lli * = l,oo, (A4a) 
||^S||ui< p|U|B|Ui. (A4b) 

In particular, ||AB||i < ||^4||oo ll-^ll i- 

Another useful fact is that the partial trace is contractive, i.e., 

UTrspOHijTr^pOll! < ||-XT|| X , (A5) 

for any operator X acting on the Hilbert space Ha <8> T~Lb- Actually this result extends to other unitarily invariant norms, with a 
prefactor depending on the dimension of the traced-out Hilbert space li47ll . 



Appendix B: Markov Approximation Bound 

We wish to derive an upper bound associated with the error from the approximation made in Eq. ( ff3| l, which involves the 
replacement of pg(t — t) by ps(t) and the extension of the upper integration limit to infinity, i.e., 

X:= / dT{A p (t-T)p s (t-T)A a (t) + ...}B a p(T) k> / dT{A fi {t-r)p s (t)A a {t) + ...}B a p(T) 
Jo Jo 

= 2 + A 1 + A 2 , (Bl) 
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where 



Ai := / dT{A fi (t-T)\p s (t)-ps(t-T)]A a (t) + ...}B a p(T) 



dr {A p (t - r)p s {t -r)A a (t) + ...} B a0 {r) , 



(B2) 



and where the ellipsis denotes the three other summands appearing in Eq. ( fT3j ). The Markov approximation is more accurate the 
smaller the error terms Ai and A2. 

Consider first the Ai term. We shall show that it is of order <? 2 tJ>. For simplicity we consider only the first of its four 
summands (the bounds for the other three are identical to that for the first): 



Jo 



d T Ap{t - r) \p s {t) -p s (t- t)] A a {t)B a p{T) 



< / dr\\Ap{t - r)|| 00 r max 

t'£[t-T,t] 



dps(t') 
dt> 



|^a(*)||oo|Ba/s(T) 



dr t max 

t'e[t-T,t] 



dpsit 1 ) 



df 



\B a0 {r)\ , 



(B3) 



where we used the triangle inequality, 1 1 ) 1 1 00 = \\A\\ = 1 and < || JTjloo ||yj|i (see Appendix [A]>. We can now use 

Eq. (|9) to upper-bound the time derivative: 



d _ , x 



<5 2 E f dT (WMt ~ r)\U\p s (t - rJUA^l 
= 4 3 2 ]T f dr\B aP (r)\ , 

. o Jo 



\B aP (T) 



(B4) 



a, fS ' 

where the factor of 4 is due to the same number of summands appearing in Eq. d9]l. Substituting this bound back into Eq. ( fB3 
we have 



dr t max 

t'e[t-T,t] 



dps(t') 



df 



a/3 



\B a p(T)\<4g 2 y2 dTT\B a p(r)\ max / dr'\B a p{T')\ 
^— ' t'e[t-T,t]J 



<ig 2 J2 dr t \B a0 (r)\T B < Ag'T^Y, 1 > < B5 > 

a/3 J° a/3 



where we used f Q di J \B a p{T')\ < f n °° dT l \B c ,s{T , )\ andEq. (H). 

This is to be compared to the term we use after the Markov approximation: 



drAp(t - T)pp s (t)A a (t)B a p(r) 



< / dr\B aP (T)\ ~t b 
1 Jo 



(B6) 



Comparing Eqs. ( |B5| l and ( |B6[ ), we see that the relative error is 0[(grs) 2 ]. 
Next consider the A2 term. We have 



drA p {t - T)p s (t - T)A a {t)B aP {T) 



< 



dr\\Ap{t - r)|U|p 5 (* - T-JIKH^WHoolB^WI 

/oo 
dr\B a p(r)\. 



(B7) 



The last integral converges provided |B a ^(r)| ~ 1/t x with x > 1. This is typically the case. E.g., for the Ohmic spin-boson 
model discussed in subsection VB we show that for t > r tr ~ (3\n(f3oj c ), the bath correlation function |B a; g(T)| ~ ffrf?- I n 
this case, then, we can set t — r tr and bound 



|A 2 ||i< 



mi 



2 eoa 



dT 



1 



V9 



(B8) 



Pu c J Tti ' r 2 /3 2 w c In()9w c ) ' 
which tends to zero as the cutoff tends to infinity at a fixed finite temperature (even if the lower integration limit is replaced by a 
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constant), as expected in the weak coupling limit assumed in our work. Note that the infinite temperature limit, where Eq. ( |B8| ) 
diverges, is incompatible with weak coupling, and requires the so-called singular coupling limit [22, 30 1. 



Appendix C: Properties of the spectral-density matrix F a p(uj) 

Introducing the Fourier transform pair 

/oo J rOO 
-oo J — oo 

and using the property that 

j~dre^> = ^-J) + iv(^) , 
where V denotes the Cauchy principal value^jwe have, using Eq. (14) , 



x r°° rli >' 

are / -^e J a p{u ) 



2tt 



in agreement with Eq. (15) , where 



in agreement with Eq. ( |16) . 
Note that: 



S a p{ui) 



r i 



-no 2?T U-W' 



B^(t,0) = (B a ( T )Bp(0))* = Tr[B a (T)B0(O)p B }* = T±[Bp(Q)B a (T)p B }=Bp a (0 > T) 
When ps commutes with H B (which we have assumed), we further have 



Bp a (0,T) = Tr [B p (0)U B (r)B a (0)U B (T)p B 
Thus the spectral-density matrix satisfies: 



= Tr 



PBU B (T)B (O)U B (T)B a (O)] = B Pa {-T,0) 



(CI) 



(C2) 



(C3) 



(C4) 



T* ap (cj) = / dr e-^Blpir, 0) = / dr e~ MT B p a (-r, 0) 
Jo Jo 



dr e iUT Bf} a (T,0) = -~/p a (uj) - iSp a (uj) 



(C5) 



(C6) 



(C7) 



Appendix D: Proof of the KMS condition 

The proof of the time-domain version of the KMS condition, Eq. (17) , is the following calculation: 

(B a { T )B b (0)) = Tr[p B U B (T,0)B a U B (T,0)B b } = ^Tv[B b e-<-^ HB B a e- lTHs ] (Dla) 

= ^Tr[B b e^+^ HB B a e^ T+ ^ HB e-^ HB ] =Tr[p B B b U B ( T + i(3,0)B a U B (T + ip,0)} (Dlb) 
= (B b {Q)B a { T + ip)) . (Die) 
Note that using the same technique it also follows that 

(B a (r)B b (0)) = (B b (-r - ip)B a (0)) . (D2) 



By definition, V ( i ) [/] = lim,,-^ / R w_ e e i dx, where / belongs to the set of smooth functions with compact support on the real line 1 
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To prove that this implies the frequency domain condition, Eq. ( p~8[ >, we compute the Fourier transform: 

/CO />oo 
dre^ T (B a {T)B b {Q)) = dre^ T (B b {-r - i0)B a {Q)) (D3) 
-CO J — CO 

To perform this integral we replace it with a contour integral in the complex plane, § c dre luJT (B b (—r — i/3)B a (0)), with the 




FIG. 9. Contour used in proof of the KMS condition. 

contour C as shown in Fig. [9] This contour integral vanishes by the Cauchy-Goursat theorem |48| since the closed contour 
encloses no poles (the correlation function (B b (T)B a (0)) is analytic in the open strip (0, —iff) and is continuous at the boundary 
of the strip [49 1), so that 



(...)=0= / (. 
c Jt 



(...)+ / (...)+ / (...) 



(D4) 



where (...) is the integrand of Eq. ( |D3| l, and the integral J is the same as in Eq. {D3}. After making the variable transformation 
t = —x — ij3, where x is real, we have 



(...) 



-e** / e-^ x (B b (x)B a (0)) 



(D5) 



Assuming that (_B a (±oo — i/3)B b (0)) = (i.e., the correlation function vanishes at infinite time), we further have /*(•••) = 
/.(•••) =0, and hence we find the result: 



dre iur {B b {-r - ip)B a {0)) 



which, together with Eq. (|D3[), proves Eq. (fT8 



e- iUT {B b (T)B a (0)) = eP»j ba (-u>) , 



(D6) 



Appendix E: Non-Adiabatic Corrections 

For completeness we provide a brief review of pertinent aspects of adiabatic perturbation theory, and the derivation of the 
adiabatic condition relating the total evolution time to the ground state gap. See, e.g., Ref. [5 1 for additional details. 



1. Adiabatic perturbation theory 

To consider adiabatic corrections, we recall that Us satisfies 

id t U s {t,t') = H s (t)U s (t,t') , 

where 

H s (t)=Y / Za(t)\E a (t))(e a (t)\. 



(El) 



(E2) 
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Define the "adiabatic intertwiner" W(t, t'): 

W(t,t') = V|e o (t))(e o (0|=T+expH / dr K(t)] , (E3) 

where the "intertwiner Hamiltonian" is 

K(t) = i[d t W{t,t')}W\t,t') =iY J \e a {t)){s a {t)\ (E4a) 

a 

= i[P {t),P (t)}, (E4b) 

and where Pa(t) = \eo(t))(so(t)\ is the projection onto the ground state of Hs(t). The result ( |E4b[ ) is by no means obvious and 
is proven in subsection |E 3| 



To extract the geometric phase we define 

H G {t) = J2Mt)\^a(t))(e a (t)\ , (E5a) 

a 

<t>a{t) = i(e a (t)\i a (t)) , (E5b) 
H' s (t) = H s (t) - H G {t) . (E5c) 

Now define V via the "adiabatic interaction picture" transformation: 

V(t,t') = W^(t,t')U s (t,t'), (E6) 

along with 

H s (t,t') = W\t,t')H s (t)W(t,t') = J2 s a(t)\s a (t'))(e a (t')\ , (E7a) 

a 

H G (t,t') = W\tX)H G {t)W(t,t') = J2<l>a(t)\za(t')){e a (lS)\ , (E7b) 

a 

H' s (t,t') = H s (t,t') - H G (t,t') , (E7c) 
K(t,t') = W\t,t')K(t)W(t,t') = iW\t,t')d t W{t,t') . (E7d) 

Note that the time dependence of Hs and H' s (t) is entirely in the energy eigenvalue and not in the eigenstates. Then V obeys 
the Schrodinger equation: 

id t V{t,t') = Hf(t,t')V{tj) . (E8) 

where 

Hf{t,t') = H s (t,t') - K(t,t') = W^(t,t')[H s (t) - K{t)]W{t,t') . (E9) 
When the evolution is nearly adiabatic H^(t, t') is a perturbation, so that we consider a solution of Eq. ( |E8| l for V of the form: 

V(t,t') = V (t,t')(t + V 1 {t,t') + ...) , (E10) 
with the zeroth order solution associated with the purely adiabatic evolution, including the geometric phase: 



V (t,t')=T + exp 



i / dTH' s (r,t') 



(Ella) 



Uf{t,t') = W(t,t')V (t,t') = £ \e a (t))(e a (t')\e- l ft' dT ^)-<p a (T)] (Enb) 

a 

Differentiating with respect to t, we have 

Ufit, t 1 ) = W(t, t')V (t, t') + W(t, t')V (t, t 1 ) (E12a) 

= -iK(t)W(t, t')V (t, t') - iW(t, t')H' s (t, t')V (t, t') (E12b) 

= -i [Hf(t) - H G (t)] Uf(t, t 1 ) , (E12c) 
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where 

Hf{t)^K{t)+H s {t). 
Plugging the Eq. ( |E10| > expansion into Eq. ( |E8| l, we obtain, to first order in 

idtV 1 (t,t') = -Vj(t,t') \k(t,t') - H G {t,t')\ V (t,t') 

= -iUf\tX)dtW{tX)V a {t,t') + uf{t,t')H G (t)Uf{t,t') ■ 
Note that Uf\t, t')H G {t)Uf{t, t') = <Pa(t)\Sa(t'))(e a {t')\. Therefore, integrating, and using 

Vba(t,t') = J druj ba (T) , uj ba (r) = [e b (r) - (j> b (t)] - [e a (r) - 4> a {t)] 



we can write the solution as: 

Vi{t,t') 



dr 



(El 3) 

(E14a) 
(E14b) 

(E15) 



uf\r, t')d T W(r, t')V (r, t>) - ]T Q (r)| £a (t')> <^«C*0l 

Therefore, the system evolution operator can be written as: 

U s (t, t') = Uf(t, t 1 ) + Q(t, t')Uf(t, t') + O (ij 2 ) , 
where the correction term can be made appropriately dimensionless in a more careful second order analysis, and where 

Q(M') = Uf{t t l/)V 1 {t,l/)Uf{t,l/)=^2e- i ^'^ (- f dre^^(e a {r)\e b (T))] \e a (t))(e b (t)\ . 

a^b V Jt ' ' 



(E16) 



(E17) 



(E18) 



2. Adiabatic timescale analysis 

Equation ( |E17| i shows that the first order correction to the purely adiabatic evolution is given by Q. Thus the adiabatic 
timescale is found by ensuring that the matrix elements of Q are all small. Let us show that a sufficient condition for this is 
-C 1 [Eq. (f27|>]. More rigorous analyses replace the <C condition with an inequality relating the same parameters to the 
fidelity between the ground state of Hg(tf) and the solution of the Schrodinger at tf, e.g., ||6l [10UT21 . 

Starting from Eq. ( |E18[ ) and taking matrix elements we have 

Q ab (t,t') = e -W*>*') (- j\Te^ T ^(e a (T)\d T \e b (T))^ . (E19) 
Changing variables to dimensionless time s = t/tf, ^ ba (t, t') becomes tf j*,^ ds' io ba {s') — tffj, ba (s, t' /tf), and 

f dTe i ^ T ^{e a {r)\d T \e b {r))= f ds' j*^'*' ^ (e a (s')\d s ,\e b (s')) , (E20) 
Jt' Jt'/tf 

where n now involves a dimensionless time integration. Using the fact that e lt/ ' ibo ( s >* l ts ' = -. — l —i—f\ -^ 7 e ltl ^ ba ( s >' l tf \ we 

1 00 t fu a b{s / ) as' 



3 To see that this is an expansion in powers of 1/t t , transform to the dimensionless time variable s = t/tf. 
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can integrate Eq. ( |E20| i by parts as 

ds > e 4t /' i -^' t '^)(e a ( s ')|a s ,|e 6 ( S ')) = 



t'/t f 



t f L0 ab (s') 
i 



*i*fM''f/*f)(e a (tr)\d a ,\E b (</)) 



t'/tf 



it f ji ab (s',t'/tf) j 

, , ds> — — JL{s a {s')\d s ,\E b {s')) 

tf Jt'/tf UabW) as' 



(E21) 



Now note that for non-degenerate energy eigenstates, differentiation with respect to t of Hs(t) \£ a (t)) — £a(t) \e a (t)), and 
substitution of s = t/tf in Eq. (j26j, directly yields the relation^] 



(e b (t)\e a (t)) = 



(e b (t)\H s (t)\e a (t)} 



U a b{t) 



Atf 



(E22) 



Substitution into Eq. ( |E21[ ) yields, with the help of Eq. ( |E22| > 

\{e a (s)\d s H s (s)\e b (s)}\ Ws 



\Qab(t,t')\ < 



U 2 a b (s)tf 



t'/tf 



t/tf pit f Labis', t'/tf) J 

ds' 2 , ^{e a { S ')\d s ,H s ( s ')\s h ( s ')) 

t'/t f ^ab^Ff ds 



Continued integration by parts will yield higher powers of the dimensionless quantity ( £ "( a )ldT#sOO|s b (s)) ^ f; nus a su ffj c i en t 
condition for the smallness of \Q a b(t, f)\ for all a, b is that 

max s;Qib \(e a (s)\d s Hs(s)\e b (s))\ 



un S ;a,bul b (s)tj 



< 1 



(E23) 



namely the condition in Eq. ( [27] ), where we assumed that the minimal Bohr frequency is the ground state gap, u>ifi. Of course, 
this argument is by no means rigorous, and in fact it has been the subject of considerable discussion, e.g., [ 50 - 55 1. For rigorous 
analyses see, e.g., Refs. 15TI71 [T0HT21 . For our purposes Eq. ( [27] ) suffices. 



3. Intertwiner Hamiltonian 

Here we provide an explicit proof of Eq. ( |E4b[ ), a well-known result due to Avron et al. [56 1. The original proof lacks some 
detail, so our purpose here is to fill in the gaps, for completeness of our presentation (see also Ref. Il25ll for details of the proof). 
Define the ground state and orthogonal projectors Po(t) = \eo(t)) (eo(t) \ and Qo(t) = I — Po{t). Recalling Eq. ( |E17| i, we have 

P (t)U s (t,t')P (t') ^Pom^it^t^Poit'j + Po^Qit^^^t^Poit^ + O^ . (E24) 

Using Eq. pTb] i we find U a s d (t,t')P Q (t') = P (*)e~*-# dTE °( T \ so that up to a phase P (t)Q(t,t')Uf(t,t')P (t') = 
Po(t)Q(t,t')Po(t'). However, Q(t,t') [Eq. ( |E18[ )] contains only off-diagonal terms [since we subtracted the Berry connec- 
tion in Eq. ( |E5c| i], so that Po(t)Q(t,t')Po(t') = 0. Thus, to first order in 1/tf the evolution generated by Hs(t) keeps the 
ground state decoupled from the excited states, i.e., 

P {t)U s {tX)P^t') = Uf{t,t')P () {t') + o{t- J 2 ) , (E25a) 

Qo{t)U s (t,t')Q (t') = Ug{t,t')Q Q {t') + o(tj 2 ^ . (E25b) 

Letting r = t — t' > and expanding around t, we then have, using Eqs. (JETJ and < |E3) , 

P (t)[l + iTHs(t)][P (t) - TP (t)] = [I + iTHt(t)][P (t) - TP (t)] + O (tf^j , (E26a) 

Q (t)[l + iTH s (t)][Q (t)-TQ (t)} = [I + irHf{t)][Qo{t)-rQ Q {t)]+o{tJ 2 ^ . (E26b) 



4 Differentiating (where x = d t x), we have Hs \e a (t))+Hs(t) \e a (t)) = e a (t ) \£ a {t))+£ a {t) \e a (t)). Taking matrix elements gives {e b (t)\ Hs \e a {t)) + 
£b{t){£b(t)\£a(t)) = £a{t)& a b + £a (t) (sb (t) \ia (t)), which yields Eq. \E22) provided |£(,(t)} is an eigenstate non-degenerate with [e a (t)), 
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Using the projector properties P = Pq => PqPo + PoPo = Pq, and [Hs(t),Po(t)] = (similarly for Q ), this simplifies to 

- rP Q (t)P (t) + iTH s (t)P (t) = -rP (t) + irHf(t)P (t) + O (^ 2 ) , (E27a) 
-tQ (*)Oo(*) + iTH s {t)Q (t) = -rQ (t) + %THf(t)Q (t) + (tf} , (E27b) 

and, after dropping the O ^/ 2 ) correct i° ns > t0 

i[Hf(t) - H s (t)]P (t) + Po(t)P Q (t) - , (E28a) 
i[Hf(t) - H 8 (t)]Q (t) + Qo(t)Q (t) = . (E28b) 

Inserting Qo(t) = I — Po(t) and Qo{t) = —Po(t) into Eq. ( |E28b[ ), and subtracting it from Eq. ( |E28a) , we find, using PoPq + 
PqPq = Po once more, the desired result: 

Hf(t) = H s (t) + i[P (t), P (t)} , (E29) 
which, together with Eq. ( |E13| l, proves Eq. ( |E4b| i. 



Appendix F: Short Time Bound 



We wish to bound the error associated with neglecting in Eq. ( |40| i, i.e., we wish to bound 

||e(i,r)|U = \\U s (t-T,0)-e iTH ^Uf(t,0)\U = \\Uf(t, 0)e"^Wt4(i, t - r)U s (t, 0) 
Using that the operator U(t) = e _iTHs(t) f/^(t, t — r) satisfies: 



dr 



U = -i[H s (t) - H s {t - t)]U(t) 



we can write a formal solution for U as: 



U(t) =1-i dt' [H s (t) - H s {t - 0] U{t') . 
Jo 



Therefore we can bound: 

l|e(*,r)iu 



i I dt'Uf(t,0) [H s (t) - H s (t - t')] U{t')U s (t,0) - I 
o 



<min{2, ||Q(* > 0)|| oo + / dt'\\[H s (t)-H s (t-t')] 
Jo 



o(t- 2 )} 



< min|2, — h - 

1 ' A 2 t f 2 t'e[t-r.t 



t 2 max WrHsit')^ + 0(tf)} 



(Fl) 

(F2) 

(F3) 

(F4a) 
(F4b) 
(F4c) 



where we used Eq. (JETT) and the fact that supoperator norm between two unitaries is always upper bounded by 2 in the second 
line, and the standard adiabatic estimate to bound \\Q(t, 0)\\oo (recall subsection E2 1. While h of Eq. |26) is expressed in terms 
of a matrix element, a more careful analysis (e.g., Ref. 1 10|) would replace this with an operator norm. Thus we shall make the 



(F5) 



plausible assumption that h ~ tf max t 'e[t~T,t] \\df Hs(t')\\oc, and, dropping the subdominant 0(tf ), we can write 



||e(t,r)||oo <min{2, 



T 2 h, 



AH, 



26 



We can now bound the error term in Eq. ( g2b) , Let X(t, r) = e lrHs ^Uf(t, 0), so that C/ 5 (t - r, 0) = r) + 9(t, r). We 
can then write 



dr^^-rJpsW^aWBa/jW = / dTX\t,T)ApX{t,T)~ Ps {t)A a {t)B a0 (T) (F6a) 

Jo 

/>oo />oo 

+ / drXt(t 1 r)^e(t,r)p s (t)A»(t)B 0i8 (T)+ / drXet(t,r)A^X(t,r)p s (t)A Q (t)B^(r) (F6b) 
Jo Jo 

+ dTQHt,T)ApQ(t,T)p s (t)A a (t)13 a0 (T), (F6c) 



The first term on the RHS of ( |F6a[ ) is the approximation we have used in Eq. ( |42b[ ). The terms in |F6b) and ( |F6c| > can be bounded 
as follows, using Eq. ( |F5[ >, the unitarity of X, the fact that ||/3si|oo < 1, and recalling that ||A a ||oo = 1. First we assume that 
Eq. ( fl2] > applies. Then: 



< 



drX* (t, r)A e(t, T)p s (t)A a {t)B a p (r) 
rfr||e(t,T)|| 0o |B^(r)|<mm{2T S , 
dTei(t,T)ApO{t,T)p s {t)A a (t)B a p{T) 



dTX&(t,T)A fj X(t,T)p S {t)A a {t)B a p{T) 



A 2 */ + t, 1 



(F7a) 



< / dr||e(t,T)||2 |B a/S (r)|<niin{47- B ,2 
Jo 



>■ (F7b) 



where in the last inequality we used the fact that if x < 2 then [min(2, x)] 2 — x 2 < 2x, and if x > 2 then again [min(2, a;)] 2 = 
2min(2,a;) < 2x, with x = ^fe - + x^' m or d er to avoid having to extend Eq. ( p~2] > to higher values of n. In all, then, the 

3 , 

approximation error in Eq. ( |42b[ > is 0[min{rs, ^fp- + "^-}]- 

Next we recall from the discussion in subsection II C that Eq. ( |T2] > must, in the case of a Markovian bath with a finite cutoff, 
be replaced by the weaker condition ( [23^ , reflecting fast decay up to r tr , followed by power-law decay. In this case the terms in 
( |F6b| i can instead be bounded as follows: 



dTX\t,T)ApQ{t,T)p s (t)A a (t)B aP {T) 



dTXQi{t,T)A X(t,T)ps(t)A a {t)B a p(r) 



< 



dT\\Q(t,T)\UB aP (T)\+ / dT\\Q(t,T)\UB a0 (r)\ 



r B h r%h. 



<min { 2. B ,^ + ^ } + 2^ 



(F8) 



in place of ( |F7b[ ). A similar modification can be computed for the term in ( |F6c[ ). To compute the order of we recall that 

2 

r tr ~ /31n(/3cj c ) ^ /3 [Eq. (|77)], and that tm = \f2p/u) c . Thus ^r 1 ~ [w e ln(/3w c )] _1 . It follows that we can safely ignore the 

2 

term in Eq. ( |F8| l provided Eq. f79| ) is satisfied. The analysis of term in ( |F6c) > does not change this conclusion. 



Appendix G: Derivation of the Schrodinger picture adiabatic master equation in Lindblad form 

Starting from Eq. (|53) and performing a transformation back to the Schrodinger picture, along with a double-sided adiabatic 
approximation, yields 

U s (t,0) f drA p {t - r)p s (t)A a (t)B a p(T)ul(t,0) + h.c. w (Gla) 
Jo 

Ap^Aa^ILu^psffill^^Tapiu) + h.c. = L uAt)ps(t)Lt )a (t)r a p(w) + h.c. . (Gib) 
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Similarly, the other terms yield 

-£/<?(*, 0) f dTA a {t)A {t-T)p s {t)B a p{T)ul{t,Q)n-Y^Ll t ^ (G2) 

Using Eqs. ( fl"5j ) and ( p"6| ) for the spectral-density matrix and its complex conjugation, we are able to combine the Hermitian 
conjugate terms, starting from the terms in Eq. ( |Glb[ ): 



Y, r Qi s(w) W*)ps(t)4,«(*) + h - 



a/3 



]T [ra/j(w)L w , / j(<)w(«)L|,, a (t) + r* /J (w)L tl ,, a (t)ps(t)4,/»(*) 



(G3a) 
(G3b) 
(G3c) 



a/3 

where in the last equality we used the freedom to interchange a and j3 under the summation sign. Likewise, the terms in Eq. ( |G2 
yield 



Yr*p(u)Ll a (t)L w ,p(t)ps(t) + h.c. 

a^ 

= £ \r aP (w)Lt ia (t)L Utfi (t)p 8 (t) +Kp(u)Ps(t)Lt tfi (t)L u , a (t) 

aj3 

= H J {4,a (*) («) , PS (*) } + *5 aj 9 M [4 , Q (t) L u ,p (t) , p S (t)] , 



(G4a) 
(G4b) 
(G4c) 



a/3 



where {, } denotes the anticommutator. Combining these results with a similar calculation for the aa and bb terms in Eqs. ( |Glb| ) 
and (|G2|i, finally results in Eqs. (|55]l and d56|. 



Appendix H: Calculations for the Spin-Boson Model 

We recall the following basic facts about the bosonic operators appearing in Eq. ( |65] l, where {n} = {n k } k is the set of 
occupation numbers of all modes: 

= 8kk', h\{n}) = ^Jnk\{n k - l,n}) , b\\{n}) = ^n k + l\{n k + l,n}) , (Hla) 
H B \{n}) = E {n} \{n}) = (£ fe n fcWfe ) \{n}) . (Hlb) 

Recalling that the bath operators fij = J2 k 9k y>\ + &fcj> we proceed to calculate (Bi(t)Bj(0)) = (Bj(t)Bj(0)) assuming 

that the bath is in a thermal Gibbs state ps = exp (—{3Hb) jZ, where Z = Tr(exp (—/3Hb)) is the partition function. We 
begin by writing: 

(BitfiBjiO)) = Tr B {u B {t^)B\U B {t^)B jP ^ 

= E (Wl4(^0)Bl|{n})({n}|C/ B (t,0)i? 3 |M)(M|p b |{™}) . (H2) 

{m},{n},{p} 

The time evolution operator acting on the eigenstates simply produces a phase, so we focus on the operator Bj 's role. 

(WI-BilW) = E (ffeVW+15{„},{ Pfc+ i,p} + fli\/pfc£{n},{ Pk -i,p}) . (H3) 

k 



bk,bl 
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Plugging this result in, we find: 



(*)*i(0)} = i £ Y. AE{mi ' E{ni)t ^ E{m} 

{m},{n},{p} k,k' 

X (^9k9k'^ n k + lyfe + l ( 5{m},{ri,n fc +l}<5{n},{p,p fc ,+l}^{p}{m} 



+ 9k9k'V n k\/Pk' + l${?n},{n,n k -l}d{n},{p,p k ,+l}5{p}{m} 
+ 9k9k'^ n k + l\ / ft'<5{m),{'',''t+l)V},(PA'-l|*W{"') 

+ 9k9i'V n kVPk'fi{7n},{n,n k ~l}(>{n},{p,p k ,~l}${p}{ m }j ■ (H4) 

The only terms that are non-zero are the middle two, giving us: 



{m} fc 



(BiffiBjiO)) =^E S f + ( e -^«»e i fa"»- J Wi>) t + e -PE lmk+1}el (E {mk+1} -E {mi )t^ (H5) 
Using that: 

OO OO - 

£ {m} - E{mk+1} =-u k , z = n e - n I— 



F —ftuJ k ' 

fc=l mfe— fc— 1 



we can write: 



(H6) 



£ ( ™ fe + l)e-^l-I%^ = r -^. (H7) 

{m} 

We can simplify our expression to 

(B^B^O)) = 1 _ ]_^ k {' '-"9l!li + e^-^gU) . (H8) 

k 

We can replace the sum with a integral: 

T2 = y2 duf(uj)S(u} - J) = / duj,J(uj) , 
i- ,.,/ ^0 Jo 



(H9) 



where /(w) is a measure of the number of oscillators at frequency ui and J(w) is the bath spectral function. Our final result is 
then: 

(A(t)^(O)) = ^T^L (e-^ff^i + e iwi "^3^i) ■ (H10) 
■/o 1 — e 

where we have assumed that oscillators with the same uj^ value interact with the i lh spin with the same interaction strength, i.e. 

g\. = g\, if w fe = u k > . (Hll) 

Plugging our result into Eq. ( [To} , we find: 

-WW" = T:'^mi (ff-..(«.l<,.<.)i e ("»-( f » + «"'"~- ( ' >l 9U.(«)l<.(,) | 8<-"»«(*))) <H12a) 

S « ( ^ (,)) =r^ra(^ P (=^) + " e """ P (=S^)) ' ,H12b) 
where denotes the Heaviside step function. 
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Appendix I: Derivation of the Ohmic Bath Correlation Function 

Here we derive the bath correlation function for the Ohmic oscillator bath with a finite frequency cutoff, Eq. ( |70| ). We start 
from the simplified expression for the spectral-density, 

ue -\u\/ul c 

TM = W 1 _ e _ /3M , (ii) 



and compute the bath correlation function by inverse Fourier transform of Eq. ( |16a| >: 

where we changed variables to x — pcj. 

The Polygamma function is defined as ip( m '> (z) = hxF(z), where T(z) — e~ x x z ~ 1 dx is the Gamma function. The 

Polygamma function may be represented for (z) > and m > as tp^ m '> (z) = (— l) m+1 J °° ^_^_ 3! dx, so that in particular, 
for m = 1, we have 

^\z)= I ^r^dx. (13) 



We can rewrite Eq. dl2} 



na 2 ( f°° xe -x[-ir/0+l+l/(Puj c )} poo -x[ir/P+l/(Pu c y\\ 

«M = ^(/ + 1 * ~ ) 

2 

= (V* 3 Hr/)9 + 1 + V(^c)) + + V(A"c))) , (Mb) 



which is the desired result. 

Appendix J: Derivation of the Effective Rate Equations in the Thermal Phase 

Here we derive the rate equation ( |85| ). 

1. Single qubit 

For illustration, let us consider a single qubit such that the Hamiltonian is given by 

H s (t) = A(t)a x - B(t)h z a z . (Jl) 

The gap at any given time in the evolution is 

A(t) = 2V 'A{tf + B(t) 2 hl . (J2) 
Since there are only two states, there are only three j(uj) terms to calculate (assuming g^ = g): 

7W -"f 2 - 7(+A(t))=W I -^ ) 7 (-A( f ))=2V ^w - ^ 

where we have taken the limit lu c — >• oo for simplicity. Let us consider the gapped phase of our evolution where t/tj is small. 
In this phase, B{t) s» 0, so we can safely assume that the energy eigenstates remain diagonal in the a x basis, with ground state 
|— ) and excited state |+), and do not change such that: 

(±|p|±) = !«±|p|±)) = p ±± . (J4) 
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The resulting equations for the density matrix elements for small times are then 

p++(t) = |(+k z |-)| 2 ( 7 (-A(t))p__ - 7 (A(t))p++) = 7*>(A(t)) (e^ A ("p- - , (J5a) 

P— (t) = 1 - />++(*) , (J5b) 

where we have used that (+|er 2 |+) = (— |er 2 |— ) = 0, (+|c z | — ) = 1, and [Hs,ps] ~ 0. Interpreting Eq. <JJ5j as a rate 
equation, we see that 7(+A) is associated with relaxation from the higher energy state to the lower energy state, and j(—A(t)) 
is associated with the excitation from the lower energy state to the higher energy state. Note that at t = 0, we assume that the 

system is initialized in a thermal state such that p++/p = e~^ A (°) w 10~ n . Since the gap remains relatively large during 

the gapped phase, the change in the initial thermal state is minimal. 



2. N qubits 

For the iV-qubit case, we can simply generalize our arguments for the single qubit case. For sufficiently small times the 
system is diagonal in the a x basis, and the lowest lying energy states can be considered to be single spin flips in the x direction. 
Furthermore, the gap between the ground state and the first excited states is A(t) = 2A(t) [Eq. ([82]>], which is large in our case. 
We can restrict ourselves to only these low lying energy states since higher excited states will have twice the gap to the ground 
state, and so their contribution will be further suppressed at short times. 

Now, using Eq. ( [83} and the KMS condition [Eq. ( [84| >] we find a similar set of equations as in the single qubit case: 

Pii « J2lapM{O\a z \i){i\a z a \O) (- Pii + e~^° Poo ) = lu (u; l0 ) (e'^poo - Pa) (J6a) 

a/3 

poo ~ 1 - Npu , (J6b) 
where we have again assumed that the initial state is the thermal state and the gap is large. 
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